import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['font.sans-serif'] = ['STHeiti']
# 使用Steinhart-Hart方程拟合热敏电阻参数 1/T = A + B*ln(R) + C*(ln(R))**3
# 由于规格书通常执提供参数B，此处使用规格书中的R-T表拟合参数
# Steinhart-Hart拟合结果:
# A = 1.1218199445e-03
# B = 2.3582428380e-04
# C = 7.7346110381e-08
# 最大温度误差: 0.1372℃
# 均方根误差: 0.0318℃

# 参考规格书 https://item.szlcsc.com/datasheet/KNTC0603%252F10KF3950/3124166.html
data = {
    'Temperature': [-40, -39, -38, -37, -36, -35, -34, -33, -32, -31, -30, -29, -28, -27, -26, -25, -24, -23, -22, -21,
                   -20, -19, -18, -17, -16, -15, -14, -13, -12, -11, -10, -9, -8, -7, -6, -5, -4, -3, -2, -1, 0, 1, 2,
                   3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
                   30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,
                   55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79,
                   80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103,
                   104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123,
                   124, 125],
    'R_Center': [345.275, 322.791, 301.925, 282.549, 264.549, 247.816, 232.254, 217.774, 204.292, 191.735, 180.032, 
                169.120, 158.941, 149.441, 140.571, 132.284, 124.522, 117.266, 110.480, 104.130, 98.185, 92.618, 87.402,
                82.513, 77.927, 73.626, 69.588, 65.797, 62.237, 58.890, 55.744, 52.786, 50.002, 47.382, 44.916, 42.592,
                40.400, 38.333, 36.385, 34.548, 32.814, 31.179, 29.636, 28.178, 26.800, 25.497, 24.263, 23.096, 21.992,
                20.947, 19.958, 19.022, 18.135, 17.294, 16.498, 15.742, 15.025, 14.345, 13.699, 13.086, 12.504, 11.951,
                11.426, 10.926, 10.452, 10.000, 9.570, 9.162, 8.773, 8.402, 8.049, 7.713, 7.393, 7.088, 6.797, 6.520,
                6.255, 6.003, 5.762, 5.532, 5.313, 5.103, 4.903, 4.711, 4.529, 4.354, 4.187, 4.027, 3.874, 3.728,
                3.588, 3.454, 3.326, 3.203, 3.086, 2.973, 2.865, 2.761, 2.662, 2.567, 2.476, 2.388, 2.304, 2.224,
                2.146, 2.072, 2.001, 1.932, 1.866, 1.803, 1.742, 1.684, 1.628, 1.574, 1.522, 1.472, 1.424, 1.378,
                1.333, 1.290, 1.249, 1.209, 1.171, 1.134, 1.099, 1.065, 1.032, 1.000, 0.969, 0.940, 0.911, 0.884,
                0.857, 0.831, 0.807, 0.783, 0.760, 0.738, 0.716, 0.695, 0.675, 0.656, 0.637, 0.619, 0.602, 0.585,
                0.569, 0.553, 0.538, 0.523, 0.508, 0.495, 0.481, 0.468, 0.456, 0.443, 0.432, 0.420, 0.409, 0.399,
                0.388, 0.378, 0.368, 0.359, 0.350, 0.341]
}

df = pd.DataFrame(data)

# 对数据进行排序（按电阻值升序）
df_sorted = df.sort_values('R_Center').reset_index(drop=True)
R_data = df_sorted['R_Center'].values * 1000  # 转换为欧姆
T_data_C = df_sorted['Temperature'].values
T_data_K = T_data_C + 273.15  # 转换为开尔文

# 定义Steinhart-Hart方程
def steinhart_hart(R, A, B, C):
    """Steinhart-Hart方程: 1/T = A + B*ln(R) + C*(ln(R))**3"""
    lnR = np.log(R)
    return A + B * lnR + C * lnR**3

# 准备拟合数据
y_fit = 1.0 / T_data_K  # 目标值：1/T

# 使用曲线拟合来估计参数
# 注意：我们使用1/T作为y值，而不是直接的温度
p0 = [1.0, 1.0, 1.0]  # 初始猜测值

try:
    popt, pcov = curve_fit(steinhart_hart, R_data, y_fit, p0=p0, maxfev=10000)
    A, B, C = popt
    
    # 计算拟合的1/T值
    y_fitted = steinhart_hart(R_data, A, B, C)
    
    # 将拟合的1/T转换回温度（开尔文）
    T_fitted_K = 1.0 / y_fitted
    
    # 转换回摄氏度
    T_fitted_C = T_fitted_K - 273.15
    
    # 计算误差
    error = T_fitted_C - T_data_C
    max_error = np.max(np.abs(error))
    rms_error = np.sqrt(np.mean(error**2))
    
    print("Steinhart-Hart拟合结果:")
    print(f"A = {A:.10e}")
    print(f"B = {B:.10e}")
    print(f"C = {C:.10e}")
    print(f"最大温度误差: {max_error:.4f}℃")
    print(f"均方根误差: {rms_error:.4f}℃")
    
    # 生成更密集的电阻值用于绘制平滑曲线
    R_fine = np.logspace(np.log10(min(R_data)), np.log10(max(R_data)), 1000)
    T_fine_K = 1.0 / steinhart_hart(R_fine, A, B, C)
    T_fine_C = T_fine_K - 273.15
    
    # 绘制结果
    plt.figure(figsize=(12, 8))
    
    # 原始数据
    plt.plot(R_data/1000, T_data_C, 'bo', markersize=4, label='原始数据', alpha=0.7)
    
    # 拟合曲线
    plt.plot(R_fine/1000, T_fine_C, 'r-', linewidth=2, label='Steinhart-Hart拟合')
    
    plt.xlabel('阻值 (Kohm)', fontsize=12)
    plt.ylabel('温度 (℃)', fontsize=12)
    plt.title('阻值-温度特性曲线 Steinhart-Hart 拟合', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.xscale('log')  # 对数坐标，因为电阻范围很大
    
    # 添加误差信息
    plt.text(0.05, 0.95, f'Steinhart-Hart方程拟合\nA = {A:.3e}\nB = {B:.3e}\nC = {C:.3e}\n最大误差: {max_error:.4f}℃', 
             transform=plt.gca().transAxes, fontsize=10, verticalalignment='top',
             bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
    
    plt.tight_layout()
    plt.show()
    
    # 绘制误差分布图
    plt.figure(figsize=(10, 6))
    plt.plot(R_data/1000, error, 'g-', linewidth=1, label='拟合误差')
    plt.axhline(y=0, color='k', linestyle='-', alpha=0.3)
    plt.axhline(y=1, color='r', linestyle='--', alpha=0.7, label='误差上限')
    plt.axhline(y=-1, color='r', linestyle='--', alpha=0.7)
    plt.fill_between(R_data/1000, -1, 1, color='red', alpha=0.1)
    plt.xlabel('阻值 (Kohm)', fontsize=12)
    plt.ylabel('拟合误差 (℃)', fontsize=12)
    plt.title('Steinhart-Hart拟合误差分布', fontsize=14)
    plt.grid(True, alpha=0.3)
    plt.legend()
    plt.xscale('log')
    plt.tight_layout()
    plt.show()
    
    # 检查是否满足误差要求
    if max_error <= 1.0:
        print("✅ Steinhart-Hart拟合满足温度误差小于1℃的要求")
    else:
        print("❌ Steinhart-Hart拟合未满足温度误差小于1℃的要求")
        
        # 如果误差太大，尝试使用简化的Beta方程
        print("\n尝试使用简化的Beta参数方程...")
        
        # Beta方程: 1/T = 1/T0 + 1/B * ln(R/R0)
        # 选择参考点（通常选择25℃）
        T0 = 25 + 273.15  # 25℃对应的开尔文温度
        R0 = 10.0 * 1000  # 25℃时的电阻值（从数据中查找）
        
        # 从数据中找到最接近25℃的电阻值
        idx_25 = np.argmin(np.abs(T_data_C - 25))
        R0 = R_data[idx_25]
        
        # 定义Beta方程
        def beta_equation(R, B):
            return 1/T0 + 1/B * np.log(R/R0)
        
        # 拟合Beta参数
        y_beta = 1.0 / T_data_K
        popt_beta, _ = curve_fit(beta_equation, R_data, y_beta, p0=[3000])
        B_beta = popt_beta[0]
        
        # 计算Beta方程拟合结果
        y_fitted_beta = beta_equation(R_data, B_beta)
        T_fitted_beta_K = 1.0 / y_fitted_beta
        T_fitted_beta_C = T_fitted_beta_K - 273.15
        error_beta = T_fitted_beta_C - T_data_C
        max_error_beta = np.max(np.abs(error_beta))
        
        print(f"Beta参数方程拟合结果:")
        print(f"B = {B_beta:.4f} K")
        print(f"最大温度误差: {max_error_beta:.4f}℃")
        
        if max_error_beta <= 1.0:
            print("✅ Beta参数方程拟合满足温度误差小于1℃的要求")
        else:
            print("❌ Beta参数方程拟合也未满足温度误差小于1℃的要求")
            
except Exception as e:
    print(f"拟合过程中出现错误: {e}")
    print("尝试使用数值方法直接求解Steinhart-Hart系数...")
    
    # 使用数值方法直接求解
    # 构建线性方程组: 1/T = A + B*lnR + C*(lnR)^3
    lnR = np.log(R_data)
    X = np.column_stack([np.ones_like(lnR), lnR, lnR**3])
    y = 1.0 / T_data_K
    
    # 使用最小二乘法求解
    coeffs, _, _, _ = np.linalg.lstsq(X, y, rcond=None)
    A, B, C = coeffs
    
    # 计算拟合结果
    y_fitted = A + B * lnR + C * lnR**3
    T_fitted_K = 1.0 / y_fitted
    T_fitted_C = T_fitted_K - 273.15
    error = T_fitted_C - T_data_C
    max_error = np.max(np.abs(error))
    
    print("使用最小二乘法直接求解Steinhart-Hart系数:")
    print(f"A = {A:.10e}")
    print(f"B = {B:.10e}")
    print(f"C = {C:.10e}")
    print(f"最大温度误差: {max_error:.4f}℃")
    
    if max_error <= 1.0:
        print("✅ 最小二乘法求解满足温度误差小于1℃的要求")
    else:
        print("❌ 最小二乘法求解也未满足温度误差小于1℃的要求")